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ABSTRACT 

When an accretion disk falls prey to the runaway instability a large portion of its 
mass is devoured by the black hole within a few dynamical times. Despite decades of 
effort, it is still unclear under what conditions such an instability can occur. The tech- 
nically most advanced relativistic simulations to date were unable to find a clear sign 
for the onset of the instability. In this work, we present three-dimensional relativis- 
tic hydrodynamics simulations of accretion disks around black holes in dynamical 
spacetime. We focus on the configurations that are expected to be particularly prone 
to the development of this instability. We demonstrate, for the first time, that the fully 
self-consistent general relativistic evolution does indeed produce a runaway instabil- 
ity 

Key words: accretion disks, runaway instability, gamma-ray bursts, numerical rela- 
tivity 



1 INTRODUCTION 

Thick and massive relativistic accretion disks around 
black holes (BHs) are thought to for m in e x treme 
core-collapse events of massive stars (Wooslev 11993 



l Abramowicz et all Il983l; iPapaloizou & Pringiel 11984 11985 



Koiimal Il98ct IWoodward et alj 19941; iFont & Daignd 12002b 



MacFadyen & Wooslevl Il999l; ISekiguchi & Shibatal 12011 
Ott et al.l201ll;l Wooslev & Hegeil2012D , and they are a normal 



outcome for the coalescense of neutron star (NS)-NS (e.g., 
Ruffert et al]|l997l; iRosswog et alj|2003l; IShibata & Taniguchil 



20061; lOechslin & Tankal |2006|; iBaiotti et all l2008t iLiu etal 
20081; iKiuchi et alj|2009t) and NS-BH binaries (e.g.,lRosswo 



2005 



2009 



Ury 

Chawla et alj|2010l;lRuffert & lanka 



Shib ata & Uryull2006t IShibata et al 



2009; 



2010 



Etienne et all 



Foucart et all 



201 ll) . Such systems are thought to be cand idates for the 
central engines of gamma-ray bursts (GRBs, Popham et al.l 
19991; IWoosievlll993l; |piranll2004 iLee & Ramirez-Ruizl |2007|; 



Mesz aros & Gehrelsl2012h 



Previous studies of the stability of accretion disks 
have shown that they can be subject to various types 
of global instabilities in a number of scenarios (e.g., 
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Zanotti et alj 120031; iTaylor et alj l201ll; iKorobkin et all 12011 
Kiuchi et al.l201ll) . Instabilitie s can result in strongly var iable 
and unstable accretion rates. lAbramowicz et al r dl983h dis- 
covered the so-called dynamical runaway instability (RI) in 
thick, self-gravitating accretion disks around BHs. The RI is 
similar to the dynamical instability in close binary systems 
that occurs when the more massive binary member over- 
flows its Roche lobe. In this case, the size of the Roche lobe 
decreases faster than the size of the binary companion, which 
can ultimately lead to the tidal disruption of the compan- 
ion and the merger of the binary system. In disks around 
BHs, a toroidal surface analogous to the Roche lobe can be 
identified. A meridional cut of this surface exhibits a cusp 
located at the Lj Lagrange point. If the disk overflows this 
toroidal Roche surface, then the mass-transfer through the 
cusp will push the cusp outwards, making a larger fraction of 
the disk matter unstable to accretion. This drives the cusp out 
even further, and leads to an exponential growth of the mass- 
transfer rate. As a result, most of the disk gets consumed by 
the BH within just a few dynamical times. 

lAbramowicz et ail fl983) found that the development 
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of the RI depends on a wide range of parameters, such 
as the disk-to-BH mass ratio Mq/Mbh an d the location of 
the inner edge of the However, their investigation 
was based on several approximations and simplifications: 
they used a polytropic equation of state (EOS) for the disk 
materi al, a pseudo-Newtonian potential to model gravity 
of BH jPaczynsky & Wiital Il980h , a disk with constant spe- 
cific angular momentum, and an approximate treatment of 
the self-gravity of the disk. Subsequent works with more 
refined approximations fo und indications of a stabilizing 
effec t due to BH rotation iWilsonl Il984l: lAbramowicz et al.1 
Il998f) , while a positive radial gradient of specific angu- 
lar mome ntum was suggested to strongly suppress the 
instability ( Abramowicz et al. 1998; Daigne & Mochkovitch 
ll997l;lFont & Daignell2002al; iDaigne & Fontll2004h . Moreover, 
studies using a Newtonian pseudopo t ential for the grav- 
ity o f the BH iKhanna & Chakrabartil Il992l; iMasuda et ail 
Il998h a nd relativistic calculations in a fixed space time back- 
ground jNishida et al.ll996l;lFont & Daignd2002bh found that 
the self-gravity of the disk aggravates the instability. Never- 
theless, the distance between the inner edge of the disk and 
the location of the cusp is probably the most crucial parame- 
ter for the development of the RI: for as long as it is too large, 
the RI is unli kely to occur. 

Recently, iMontero et alj |2010j) performed the first fully 
general relativistic simulations of thick accretion disks 
around BHs in axisymmetry for a few dynamical times. For 
the particular models they studied, they found no signature 
of a RI during the simulated time. However, the inner sur- 
face of their disk models was located away from the Roche 
surfaces (R Montero, private communication), so that the in- 
stability might not have had sufficiently favourable condi- 
tions or sufficient time to develop within the time period of 
disk evolution that was considered. Therefore, their results 
do not rule out the existence of the RI in disk models where 
the inn er edge is located clos er to the cusp. 

In lKorobkin et al. I feoill) . some of us have analyzed the 
stability of slender and moderately slender accretion disks 
around BHs with M D /M BH in the range of (0.11, 0.24) using 
three-dimensional (3D) numerical simulations in full general 
relativity (GR). Although we did observe the development 
of several non-axisymmetric instabilities in our models, we 
found no traces of the RI. This result is, perhaps, not surpris- 
ing since the inner radii of our disk models were located sig- 
nifica ntly away from the Roche surface (see iKorobkin et al.l 
feoill) for more details); 

In a similar study, iKiuchi et al.l J201lh performed simu- 
lations of self-gravitating disks around BHs in full GR with 
the aim of obtaining the gravitational wave signal from the 
non-axisymmetric Papaloizou-Pringle instability. They con- 
sidered four disk models with constant and non-constant 
specific angular momentum and disk-to-BH mass ratios of 
0.06 and 0.10. No runaway instability was observed in their 
simulations, but no information was given regarding the rel- 
ative location of the inner edge of the disk and the Roche 



surface. Therefore, it is difficult to judge if their models were 
sufficiently susceptible to develop the RI. 

In order to understand whether the RI can occur at all 
in the most general case, one first has to expore whether it 
can develop in configurations that are particularly prone to the 
instability. Such systems contain disks that exactly fill their 
Roche lobes, have significant fractions of the BH masses, con- 
stant specific angular momentum profiles and non-rotating 
BH^. It is left to future studies to explore under which cir- 
cumstances such configurations would form in Nature. Our 
aim here is to explore whether the RI can occur in principle if 
the fully dynamical general-relativistic effects are taken into 
account properly. Our 3D general relativistic simulations in- 
deed confirm that the RI occurs in these cases. 

This paper is organized as follows. In Section [2] we de- 
scribe the numerical methods used in our simulations. In Sec- 
tion [3] we describe the initial disk model. In Section [4] we 
present our results, and in Section|5]present our conclusions 
for future work. Unless otherwise noted, throughout the pa- 
per we measure distances in the units of the Schwarzschild 
radius r g = 2GM BH (0)/c 2 of the black hole at t = 0, while 
the rest of the quantities are measured in cgs units. 



2 NUMERICAL METHODS 

The numerical time-evolution in our stu dy is performed 
with the GR hydrodynamics code Thor JZink et alJ|2008l) 
and the spacetime evolution code QuiLT JPazos et alj|2007l) . 
These two codes are based on and communicate with 
each other vi a the open-source Cactus computationa l 
infrastructure iGoodale et all 120031; ISchnetter et all |2007|) . 
and us e the Carpet mesh refinement and mu ltiblock 
driver JSchnetter etaTI |2004|; ISchnetter etaTI 120061) . Initial 
models of self-gravitating disks around BHs are constructed 
using an improved version of the RNS code (Stergioulas 
120111) . We use a multiblock grid with an adapted resol ution 
which is similar to the one used in lKorobkin et al.H201ll) . The 
grid resolution is A m i n ~ 0.022 r g near the BH, A max ~ 0.2 r g 
at the disk density maximum, and A ou t ~ 5.0 r g at the outer 
boundary. For more details of the numerical methods as well 
as the gr id structure used in our simulations, we refer the 
reader to lKorobkinet alj fcoill) . 



3 INITIAL DISK MODEL 

In this study, we focus on two initial disk models, denoted 
A and B. Table Q] lists the parameters of these models, while 
Fig. rjj illustrates the structure of initial model A. We con- 
struct our disks using a polytropic EOS p = Kp T with T = 
4/3, while the time-evolution is performed with the two- 
parameter ideal gas EOS (with the same T), which allows 
for shock heating. We do not include additional physics such 
as magnetic fields, nuclear or neutrino processes. Therefore, 
the results can be rescaled using one free parameter such as 
the initial BH mass. In particular, the cgs parameters listed 



1 Here, the term "disk" refers to initial equilibrium disk configu- 
rations. Therefore, the concept of the inner edge for such systems is 
well defined. Note, however, that disk tend to spread on a viscous 
timescale. 



2 BH+disk systems with slowly rotating BHs are unlikely to be dras- 
tically different from non-rotating BH cases. Nevertheless, the latter 
are expected to be more susceptible to the RI, that is why they are in 
the focus of our work. 
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Model 


A 


B 


Specific angular momentum I [Mbh] 


3.91 


3.84 


Polytropic constant K [10 14 cm 3 g _1 ^ 3 s -2 ] 


5.35 


5.72 


Maximum density p m ax [10 12 g/cm 3 ] 


1.63 


1.33 


Disk-to-BH mass ratio M D /M BH 


0.2097 


0.1628 


Kinetic to potential energy T/\U\ 


0.2666 


0.2469 


Inner radius [r g ] 


2.082 


2.127 


Location of the cusp r L [r g ] 


2.024 


2.148 


Outer radius routing] 


12.08 


12.37 


Central radius r max [r g ] 


4.687 


4.342 


Orbital frequency at r max , O max [s _1 ] 


1414 


1581 


Orbital period at r max , P max [s] 


4.44 ■ 10~ 3 


3.96 ■ KT 3 


Potential gap AW := W fa - W L 


-1.2-10- 4 


+0.01 



Table 1. Physical parameters of the initial disk models used. Notice 
that the specific angular momentum and AW are given in units of 
c = G = l. 
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Figure 2. Evolution of the disk. Top: radial location of the cusp and 
the disk density maximum, normalized to their initial values. Centre: 
evolution of the disk rest mass and angular momentum for mod- 
els A and B in full GR as well as for model A in the Cowling approx- 
imation. For clarity, the deviations from 1 of the mass and angular 
momentum in Cowling are multiplied by a factor of 10. Bottom: Nor- 
malized amplitudes of non-axisymmetric m = 1, 2, 3 deformations 
in model A and m — 1 deformation in model B. 



Figure 1. The structure of the initial disk model A. Red and green 
lines show the contours of the effective potential W = log \iit | and of 
the disk pressure, respectively. Here, lit is the covariant f-component 
of the 4-velocity of a test particle on a circular orbit with the given 
specific angular momentum for model A. The values of W are given 
in geometrized units G = c = 1. Contours of the pressure are equally 
spaced by 0.5 in decimal logarithmic scale, starting from the max- 
imum pressure (green dot on the r-axis). The blue contour corre- 
sponds to W = 0, and the black contour marks the location of the 
toroidal Roche lobe with the cusp. The disk fills the Roche lobe al- 
most up to the cusp. The gray area marks the location of the BH 
event horizon. 



in Table Q] assume A^bh = 5 Mq. The disk-to-BH mass 
ratio A^disk/^BH is ~ 0.21 for model A and « 0.16 for 
model B. Both disks have a constant specific angular mo- 
mentum profile. The black contour in Fig. Q] corresponds to 
the Roche surface, while the green ones show the contours 
of constant pressure. The disk in model A fills its Roche lobe 
almost entirely, with a remaining gap in the effective poten- 
tial between the inner edge of the disk and the cusp of only 
AW = —1.2 x 10~ 4 . By contrast, model B is constructed to 
slightly overfill its Roche lobe by the value AW = 0.01. This 
is done in order to induce the onset of the RI right from the 
very beginning of the numerical evolution. 



4 RESULTS 

Models A and B both develop the runaway instability. There- 
fore, in the following, we concentrate on the results for model 
A, while model B will be discussed later in this section. The 
top panel of Fig. [2] shows the time evolution of the cusp ra- 
dius and the radial coordinate of the location of the maxi- 
mum disk density r max . Due to initial metric perturbations 
induced by matching the initial da ta to a vacuum BH m etric 
near the horizon (cf. discussion in lKorobkin et al.ll201lh . the 
BH mass and cusp radius settle to a new, ~ 3% smaller value 
within t < 0.5t olb . The smaller gravitational pull of the BH 
leads to a rapid increase of r max by ~ 9% within the same 
time interval. Since the cusp is now located at a smaller ra- 
dius, the disk is less likely to become subject to the RI. How- 
ever, the initial metric perturbations induce oscillations in the 
disk, which are particularly evident in the evolution of r max 
at f < 6t olb . These oscillations lead to occasional crossing of 
the Roche lobe by the disk, resulting in a steady and slow 
accretion of the disk material onto the BH at t < 9i or b- This 
is visible in the evolution of the disk mass and angular mo- 
mentum shown in the center panel of Fig. [2] 

The radius r max gradually decreases due to this slow ac- 
cretion up to the time when the inner disk radius becomes 
as small as the cusp radius. This occurs at t ~ 9t olb . At 
that point, the cusp radius starts increasing, leading to an 
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Figure 3. Colourmaps of the disk density in the vertical plane at different times in evolution of model A. Here, f or b is the orbital period at the 
maximum density radius r max . 



acceleration of the accretion. By the end of our simulation 
(f = 18i or b), ~ 75% of the initial disk mass has been accreted 
onto the BH. Such an accelerated accretion due to dynamical 
migration of the cusp towards the disk during which most 
of the disk material is swallowed by the BH is exactly the 
defining property of the runaway instability. Thus, our sim- 
ulations show that the RI can indeed occur, at least for the 
most susceptible models, which are considered here. 

It is interesting to note that the development of the ax- 
isymmetric runaway instability triggers non-axisymmetric 
deformations of the disk. The bottom panel of Fig. |2] shows 
the normalized a mplitude of non-axisy mmetric m = 1, 2, 3 
deformations (see lKorobkin et all J201lh for the exact defini- 
tion of the amplitude). The m = 1 deformation grows expo- 
nentially starting at ~ 7t olb until ~ 14i or b, reaches its peak 
value of ~ 0.027, at which point it saturates. This deforma- 
tion develops due to the so-called Papaloiz ou and Pringle in- 
stability (PPI) JPapal oizou & Pringldll984l) , e nhanced by an 
eccen tric motion of the BH (as described in iKorobkin et al.l 
l20ll . The deformations corresponding to the other values 
of m do not show a strong growth and remain below ~ 10~ 3 . 
Since both the PPI and RI develop roughly at the same time 
in our simulations, they could, in principle, interact non- 
linearly, when sufficiently large amplitudes are reached. In 
particular, the PPI-induced deformations might be responsi- 



ble for the apparent saturation of the RI towards the end of 
our simulations. Such deformations redistribute angular mo- 
mentum outwards, which, in t urn, inhibits the development 
of the RI jPaigne & Fonj|20o3) (while another factor respon- 
sible for the eventual saturation of the RI is, of course, the de- 
pletion of the disk mass). Nevertheless, the rather moderate 
amplitude of the m = 1 mode of ~ 0.027 is unlikely to drasti- 
cally affect the evolution of the RI, especially during the early 
evolution. A more detailed analysis of the non-linear interac- 
tion between the PPI and the RI is required to clarify its role. 

Figure [3] shows four snapshots of the disk density in the 
meridional plane, for model A, corresponding to the times 
^ / ^orb = 0.5, 6, 12 and 18. The snapshot at t = 0.5 t or b shows 
the disk in a perturbed state after the passage of the initial 
metric perturbation. As the disk moves away from the BH, 
several surface waves can be seen propagating to the outer 
side of the disk in an axisymmetric manner. As noted above, 
such waves have sufficiently high amplitude to overfill the 
Roche lobe and support a small amount of accretion. The 
next snapshot at t = 6 t orb shows a wider accretion stream 
and an extended structure of the outer parts of the disk, 
caused by heating by the shocks formed during the radial 
disk oscillations. Such shock heating efficiently damps out 
the radial disk oscillation and causes the disk to heat up and 
expand, similarly to what was observed in IKorobkin et alj 
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feOlli) . The last two snapshots show the disk during the rapid 
accretion phase at t = 12 t orb and t = 18 t orb . By t = 18 t orb , 
the accretion stream is very wide, but the high-density parts 
have moved closer to the black hole and occupy a smaller 
volume, indicating a smaller disk mass. 

To secure our findings against the possibility of an error 
in the treatment of hydrodynamics, we have complemented 
the fully dynamical GR simulation of model A with one 
in the Cowling approximation (i.e., on a fixed metric back- 
ground). Such a model should undergo at most a slow and 
steady accretion, in a way similar to the early (t < 6f or b) evo- 
lution of the fully dynamical model. There should be no un- 
stable growth of accretion on a fixed metric background. As 
we can see from the time evolution of the disk mass and an- 
gular momentum shown in the center panel of Fig. [2] this is 
indeed the case. 

Interestingly, model A does not develop the PPI in our 
Cowling simulation, where amplitudes of non-axisymmetric 
deformations remain below ~ 1CP 5 throughout the evolu- 
tion. Since the PPI in Cowling approximation relies on the 
existence of the non -accreting inner edge of the disk jBlaesI 
ll987l;lHawley1ll99ll) . the absence of PPI is probably caused 
by the accretion at a small rate through the cusp, which can 
completely free ze the development of PPI modes in Cowling 
approximation telaeslll987l; lHawleylll99l1) . This is different 
in the case of a dynamical metric, where the PPI can b e en- 
hanced by the motion of the BH feorobkin et alj 1201 lh . For 
this reason, accretion fails to suppress the PPI in model A 
with a dynamical metric (this is also observed in Fig. 1 of 
iKiuchi et al.l j201ll) for the models with constant specific an- 
gular momentum). 

Compared to model A, in model B the disk overflows its 
Roche lobe by a small amount AW = 0.01. In this case, the ex- 
ponential accretion should commence immediately and con- 
tinue until a significant fraction of the disk is accreted onto 
the BH within a few dynamical times. This is what we ob- 
serve in our simulation, as can be seen from the evolution 
of the disk mass and angular momentum shown in the cen- 
ter panel of Fig. [2] Notice, however, that in this case the de- 
velopment of the instability happens more slowly, because 
of the lower disk-to-BH mass ratio. The amplitude of non- 
axisymmetric deformations is also lower. In particular, the 
m = 1 mode exhibits slow and irregular growth and stays 
below ~ 10~ 3 throughout the simulation (cf. the center panel 
of Fig. [2). In the absence of any significant non-axisymmetic 
deformation, the accretion of a substantial fraction of the disk 
mass onto the BH within a few dynamical times can only be 
attributed to the development of t he RI. 

As mentioned in Section [TJ iMontero et alj i2010t) pre- 
sented axisymmetric simulations of two disk models around 
BHs with a constant specific angular momentum. Their sim- 
ulations were also fully general relativistic with dynamical 
spacetime e volution, similar to ou r simulations. The absence 
of the RI in lMontero et alj teOlOl) must be attributed to dif- 
ferences in the initial disk models (compared to our config- 
urations) and specifically to the fact that their initial models 
did not exactly fill the Roche lobe. Because of this, the in- 
stability is less likely to occur within the limited simulation 
ti me. The same argum ent may explain the absence of the RI 
m IKiuchi et alj feoill) . If confirmed, this would further un- 
derline the importance of the distance between the cusp and 
the inner disk surface for the development of RI. 



Runaway Instability 5 

5 CONCLUSION 

In this study, we have, for the first time, demonstrated that 
the runaway instability in self-gravitating accretion disks 
does indeed occur in fully dynamical general relativistic evo- 
lutions. We have selected two models that are particularly 
prone to the development of such an instability. Our models 
have an appreciable disk-to-BH mass ratio (0.21 for model A 
and 0.16 for model B, see Table [1}, a constant profile of spe- 
cific angular momentum and a non-rotating BH. Moreover, 
and perhaps most importantly for the development of the 
runaway instability, our disk model A almost exactly fills its 
Roche lobe, while model B slightly overfills it. 

Our simulations show that both models develop the 
runaway instability, exhibiting unstable accretion of the disk 
matter onto the BH within just a few dynamical times. More 
than half of the disk mass is absorbed by the BH by the end 
of our simulations, which were terminated only because we 
ran out of available computing time. 

Our results demonstrate that the runaway instability 
does indeed occur, at least in the models considered here. 
Future research will need to investigate how this depends on 
the parameters of the initial disk models, such as the disk-to- 
BH mass ratio, the gradient of the specific angular momen- 
tum of the disk, and the BH spin, in order to establish the 
astrophysical significance of the instability. 
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